OHBM2024 Tutorial: TAC Data Hands-on

Setup

Packages

First we load up the packages we’ll be working with

# remotes::install_github("mathesong/bloodstream")
# remotes::install_github("mathesong/kinfitr")

library(tidyverse)
library(bloodstream)
library(kinfitr)
library(job)

theme_set(theme_light())

Data

Here is an optional chunk to download the minified dataset. However, it should already be included in the binder.

download.file("https://www.dropbox.com/scl/fi/ax1t3lfop3pp3h4cir0m7/ds004869.zip?rlkey=6z7e9wopud9ngu26ekbd3kfks&dl=1", 
              destfile = "ds004869.zip")

unzip(zipfile = "ds004869.zip")

Blood Analysis

Here we can run a basic blood analysis. For anyone who wants to run it, I’ve placed it in a job call, which runs it in the background, but it may or may not make your session slow.

job({ 
  bloodstream("ds004869/") 
})

Alternatively, if you want to use a configuration file, this is how you could do so. You can inspect the config file in the relevant directory. But best not to run this: this takes a little longer, and it is not in a job call either, so you won’t be able to do anything until it is completed.

bloodstream("ds004869/", configpath = "ds004869/code/config_tutorial.json") 

Loading Data

Loading bloodstream outputs

We want the arterial input function data, which are called “input”.

bloodstream_data <- bloodstream_import_inputfunctions("ds004869/derivatives/bloodstreamtutorial/") %>%
  select(-measurement)

Let’s plot a few of them to see. These contain the interpolated curves, and can be used with kinfitr to perform kinetic modelling.

head(bloodstream_data, n = 6)

bloodstream_data %>% 
  slice(1:6) %>% 
  pull(input) %>% 
  map(plot)

[[1]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[2]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[3]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[4]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[5]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

[[6]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).

Loading petprep outputs (TACs and morphology)

Now we want to load up the petprep outputs. There will eventually be a function for automating the next few steps, but the PET BIDS Derivatives guidelines are not yet nailed down, so we’ll do it manually for now.

First we parse the files and load in the data for the tacs and the morphology.

petprep_data <- kinfitr::bids_parse_files("ds004869/derivatives/petprep_extract_tacs/") %>% 
  unnest(filedata) %>% 
  filter(str_detect(path_absolute, "gtmseg"))

tacs <- petprep_data %>% 
  filter(measurement=="tacs") %>% 
  filter(is.na(pvc)) %>% 
  mutate(tacdata = map(path_absolute, ~read_delim(.x, delim="\t", 
                                                  show_col_types = FALSE)))

morphdata <- petprep_data %>%
  filter(measurement=="morph") %>%
  mutate(morphdata = map(path_absolute, ~read_delim(.x, delim="\t",
                                                  show_col_types = FALSE)))

Now we have lots of little regions of interest, but they are all very small. Typically, we would want to combine regions into larger regions that we are interested in. We do so using volume-weighted averaging, which requires both the morphology and the TAC information. This will also soon be automated once the PET BIDS Derivatives are more finalised.

Creating combined TACs

We’ll create a frontal cortex region, a striatal region and a hippocampus-amygdala region.

First, let’s define the regions by their constitutent parts.

frontal_regions <- morphdata$morphdata[[1]] %>% 
  filter(str_detect(name, "frontal")) %>% 
  pull(name)

frontal_regions

 [1] "ctx-lh-caudalmiddlefrontal"  "ctx-lh-lateralorbitofrontal" "ctx-lh-medialorbitofrontal"  "ctx-lh-rostralmiddlefrontal"
 [5] "ctx-lh-superiorfrontal"      "ctx-lh-frontalpole"          "ctx-rh-caudalmiddlefrontal"  "ctx-rh-lateralorbitofrontal"
 [9] "ctx-rh-medialorbitofrontal"  "ctx-rh-rostralmiddlefrontal" "ctx-rh-superiorfrontal"      "ctx-rh-frontalpole"         

striatal_regions <- morphdata$morphdata[[1]] %>% 
  filter(str_detect(name, "Putamen|Accumbens|Caudate")) %>% 
  pull(name)

striatal_regions

[1] "Left-Caudate"         "Left-Putamen"         "Left-Accumbens-area"  "Right-Caudate"        "Right-Putamen"        "Right-Accumbens-area"

hipamg_regions <- morphdata$morphdata[[1]] %>% 
  filter(str_detect(name, "Hippocampus|Amygdala")) %>% 
  pull(name)

hipamg_regions

[1] "Left-Hippocampus"  "Left-Amygdala"     "Right-Hippocampus" "Right-Amygdala"   

Now we combine the TACs and the region sizes side-by-side

selected_tacs <- select(tacs, c(ses:rec, tacdata)) %>% 
  inner_join(select(morphdata, c(ses:rec, morphdata))) %>% 
  group_by(sub, ses) %>% 
  mutate(tacdata = map(tacdata, ~pivot_longer(.x, 
                                              cols = `Left-Cerebral-White-Matter`:`ctx-rh-insula`, 
                                              names_to = "name", values_to = "TAC"))) %>% 
  mutate(tacdata = map2(tacdata, morphdata, ~inner_join(.x, .y, by="name")))

Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`

Then we perform the volume-weighted averaging itself. I define a function to do it for one region, and another function to do it for multiple regions.

volume_weighted_average_tac <- function(tacdata, regions, regionname) {
  
  tacdata_combined <- tacdata %>% 
    filter(name %in% regions) %>% 
    group_by(frame_start, frame_end) %>% 
    mutate(volume_tot = sum(`volume-mm3`),
           volume_frac = `volume-mm3` / volume_tot,
           TAC_frac = TAC * volume_frac) %>% 
    summarise(!!regionname := sum(TAC_frac), 
              .groups = "keep") %>% 
    ungroup()
  
  return(tacdata_combined)
  
}

volume_weighted_average_tacs <- function(tacdata, regions_list) {
  
  regionnames <- names(regions_list)
  
  out <- map2(regions_list, regionnames, ~volume_weighted_average_tac(tacdata, .x, .y))
  
  out <- purrr::reduce(out, dplyr::inner_join, by = c("frame_start",
                                                      "frame_end"))
  
  return(out)
}

regions_list <- list(
  Frontal = frontal_regions,
  Striatum = striatal_regions,
  HippAmg = hipamg_regions
)


selected_tacs <- selected_tacs %>% 
  group_by(sub, ses) %>% 
  mutate(selected_tacdata = map(tacdata, ~volume_weighted_average_tacs(.x, 
                                                                       regions_list)))
  

Combining the TAC and blood (input) data

Now we combine the blood data and PET data so that we can perform modelling using both of them.

modeldata <- selected_tacs %>% 
  select(ses:rec, tacdata = selected_tacdata) %>% 
  inner_join(bloodstream_data)

Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`

Fitting TACs

Preparation

Many of these steps will be automated when the PET BIDS Derivatives are finalised. But for now, we’ll do these ourselves.

Correcting Units

The TAC data are in Bq/mL, and the bloodstream data are in kBq/mL. This is because Bq/mL is hard to read, but it is the SI unit. Future functions within kinfitr will automatically convert them.

modeldata <- modeldata %>% 
  mutate(tacdata = map(tacdata, ~.x %>% 
                         mutate(across(.cols = Frontal:HippAmg, 
                                       ~unit_convert(.x, 
                                                     from_units = "Bq", 
                                                     to_units = "kBq")))))

Adding frame midtimes and durations

We only have frame_start and frame_end, but we want the midtimes and durations for modelling. We also want the times in minutes for modelling, so we’ll convert them.

modeldata <- modeldata %>% 
  mutate(tacdata = map(tacdata, ~.x %>% 
                         mutate(frame_start = frame_start / 60,
                                frame_end = frame_end / 60) %>% 
                         mutate(frame_dur = frame_end - frame_start,
                                frame_mid = frame_start + 0.5*frame_dur)))

Adding weights

There are many different weighting functions, many of which are included in kinfitr. I’ll use the default one for now. We usually run the weights on a large global grey matter brain region. For today, we’ll just take a mean TAC of our three regions and use that.

modeldata <- modeldata %>% 
  mutate(tacdata = map(tacdata, ~.x %>% 
                         mutate(meanTAC = rowMeans( .x %>% select(Frontal:HippAmg) )) %>% 
                         mutate(weights = weights_create(t_start = frame_start,
                                                         t_end = frame_end,
                                                         tac = meanTAC))))

Fitting a single TAC

Now we will fit the two-tissue compartment model to a single TAC just to see how it looks.

First we should fi

fit_2tc <- twotcm(
       t_tac = modeldata$tacdata[[1]]$frame_mid, 
       tac = modeldata$tacdata[[1]]$Frontal, 
       input = modeldata$input[[1]], 
       weights = modeldata$tacdata[[1]]$weights, 
       vB = 0.05, multstart_iter = 5)

plot(fit_2tc)

fit_2tc$par

That worked!

Notice that in the parameters, we have estimated an inpshift. This is the delay: a shift in the timing of the blood (input) data to the same time as the PET. This is necessary both because the tracer arrives at the brain at a slightly different time from when it arrives at the arm where we measure it, but also because it is hard for these two clocks to be perfectly synchronised.

We would usually fit the delay using the first 5 or 10 minutes of the acquisition using a large region for each individual with a one-tissue compartment model and check that. Then we would use that value when fitting the later models. That information is covered more fully in the full on my website.

Fitting multiple TACs

Here we’ll use a linearised model because they fit more quickly, in this case Logan. They are also less sensitive to the delay time, so we can probably reasonably safely ignore it. But most linearised models require a t* time to operate.

Selecting a t* time

Let’s choose an appropriate t* value

modeldata %>% 
  ungroup() %>% 
  filter(ses=="baseline") %>% 
  slice(1:5) %>% 
  mutate(tstarplot = map2(tacdata, input, 
     ~Logan_tstar(
         t_tac = .x$frame_mid, 
         lowroi =  .x$HippAmg,
         medroi =  .x$Frontal,
         highroi = .x$Striatum, 
         input = .y,
         vB = 0.05)
     )) %>% 
  pull(tstarplot)

Warning: There were 15 warnings in `mutate()`.
The first warning was:
ℹ In argument: `tstarplot = map2(...)`.
Caused by warning:
! Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 14 remaining warnings.

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

Ok, let’s use 13 frames, corresponding with 57.5 minutes.

Fitting

Let’s focus on the frontal cortex, and run the fitting for all measurements.

modeldata <- modeldata %>% 
  group_by(sub, ses) %>% 
  mutate(Logan_fit = map2(tacdata, input, ~Loganplot(t_tac = .x$frame_mid,
                                             tac = .x$Frontal, 
                                             input= .y, 
                                             tstarIncludedFrames = 13, 
                                             weights = .x$weights, 
                                             vB = 0.05, 
                                             dur = .x$frame_dur)))

Plotting

Let’s see a few fits

map(modeldata[1:6,]$Logan_fit, plot)

[[1]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

[[2]]

Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).

[[3]]

Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).

[[4]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

[[5]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

[[6]]

Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).

Outcomes

Now let’s observe the resulting V_T values.

Logan_outcomes <- modeldata %>% 
  select(sub, ses, Logan_fit) %>% 
  mutate(Vt = map_dbl(Logan_fit, c("par", "Vt"))) %>% 
  select(-Logan_fit)

ggplot(Logan_outcomes, aes(x=ses, y=Vt, colour=sub, group=sub)) +
  geom_point(size=3, colour="black") +
  geom_point(size=2.5) +
  geom_line() + 
  scale_y_log10()

As we can see, they are reduced after blocking with calecoxib as expected.